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Seven-transmembrane receptors (7TMRs) are involved in nearly all aspects of chemical 
communications and represent major drug targets. 7TMRs transmit their signals not only via 
heterotrimeric G proteins but also through p-arrestins, whose recruitment to the activated receptor 
is regulated by G protein-coupled receptor kinases (GRKs) . In this paper, we combined experimental 
approaches with computational modeling to decipher the molecular mechanisms as well as the 
hidden dynamics governing extracellular signal-regulated kinase (ERK) activation by the 
angiotensin II type 1A receptor (AT 1A R) in human embryonic kidney (HEK)293 cells. We built an 
abstracted ordinary differential equations (ODE) -based model that captured the available knowl- 
edge and experimental data. We inferred the unknown parameters by simultaneously fitting 
experimental data generated in both control and perturbed conditions. We demonstrate that, in 
addition to its well-established function in the desensitization of G-protein activation, GRK2 exerts a 
strong negative effect on p-arrestin-dependent signaling through its competition with GRK5 and 6 
for receptor phosphorylation. Importantly, we experimentally confirmed the validity of this novel 
GRK2-dependent mechanism in both primary vascular smooth muscle cells naturally expressing 
the AT 1A R, and HEK293 cells expressing other 7TMRs. 
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Introduction 

Seven transmembrane receptors (7TMRs) regulate nearly 
every known physiologic processes in mammals. They 
represent the largest class of cell surface receptors and are 
the target of up to 40% of current pharmaceuticals (Ma and 
Zemmel, 2002) . Signaling by 7TMRs is classically mediated by 
receptor coupling to heterotrimeric G proteins that activate a 
variety of effectors, including second messengers and the 
mitogen-activated protein kinase (MAPK) cascades (Reiter and 
Lefkowitz, 2006) . G-protein coupling is rapidly hampered by a 
two step process, which begins with phosphorylation of the 
agonist-occupied receptor by G protein-coupled receptor 
kinases (GRKs; Lefkowitz, 1998). Cytoplasmic (3-arrestins are 
subsequently recruited to the GRK-phosphorylated receptor 
and have key roles in both receptor desensitization and 
internalization (Goodman et al 1996; Lefkowitz, 1998). 

Beside this classical paradigm, a growing body of evidence 
has revealed that p-arrestins also serve as signal transducers 
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and adaptors (Lefkowitz and Shenoy, 2005; Reiter and 
Lefkowitz, 2006). So far, the best-characterized signaling 
mechanism to be stimulated by p-arrestins is the extracellular 
signal-regulated kinase (ERK) signaling pathway. p-Arrestins 
have been shown to scaffold components of this MAPK 
cascade (i.e., Raf-1, MEK1 and ERK) in complexes with 
receptors, leading to activation of ERK (DeFea et al 2000; 
Luttrell et al 2001; Ahn et al 2003). Inhibition of p-arrestin 2, 
using specific small interfering RNA (siRNA), impairs the 
angiotensin II type 1A receptor (AT 1A R) -stimulated ERK 
signaling, while p-arrestin 1 depletion enhances AT 1A R- 
mediated ERK activation (Ahn et al, 2004b). Spatial distribu- 
tion of ERK activated by the G protein- and the p-arrestin- 
dependent mechanisms is distinct: the G protein-dependent 
pathway triggers the nuclear translocation of pERK while the 
P-arrestin 2-activated ERK is confined to the cytoplasm (DeFea 
et al 2000; Luttrell et al 2001; Ahn et al 2004a). The two 
pathways also have distinct activation kinetics: the rapid and 
transient G protein-dependent activation and the slower but 
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persistent (3-arrestin 2-mediated activation (Ahn et al, 2004a). 
Interestingly, transient and sustained ERK activation have 
been shown to regulate cell fates such as growth and 
differentiation (Sasagawa et al, 2005). In addition, the 
(3-arrestin 2-dependent ERK pathway can be activated inde- 
pendently of G proteins as shown using the (DRY/AAY) 
mutant of AT 1A R or a modified angiotensin II peptide (SII) 
(Wei et al, 2003). 

Noteworthy, different GRK subtypes have been shown to 
have specialized regulatory functions for the G protein and 
(3-arrestin-dependent signaling mechanisms. Activation of the 
(3-arrestin 2-dependent ERK pathway by the AT 1A R (Kim et al, 
2005), V2 vasopressin (Ren et al, 2005), (3 2 adrenergic (Shenoy 
et al, 2006) and follicle-stimulating hormone (FSHR) (Kara 
et al, 2006) receptors specifically requires GRK5 and GRK6 
action. Second messenger generation by both V2 vasopressin 
(Ren et al, 2005) and HI histamine (Iwata et al, 2005) receptors 
is negatively regulated by GRK2 but unaffected by GRK5 or 6. 

A number of studies have provided mathematical models of 
ERK activation by EGF (Kholodenko et al, 1999; Orton et al, 
2005; Borisov et al, 2009; Kholodenko et al, 2010), including 
the dynamics of transient and sustained ERK activation 
(Sasagawa et al, 2005; Nakakuki et al, 2010). ERK response 
to a G protein-coupled 7TMR in yeast and mammals (Hao et al, 
2007; Csercsik et al, 2008) as well as the role of GRK in the 
desensitization, internalization and recycling of the (3 2 
adrenergic receptor ((32 AR) (Violin et al, 2008; Vayttaden 
et al, 2010) have also been modeled. In addition, a number of 
other studies have modeled different aspects of 7TMRs 
signaling, including calcium signaling (see Linderman, 2009 
for a recent review). However, no mathematical model has 
been developed to address GRK-mediated regulation of ERK 
activation by a 7TMR. In this study, we used experimental and 
computational modeling approaches to decipher the molecu- 
lar mechanisms governing ERK activation by the AT 1A R in 
human embryonic kidney (HEK)293 cells. 

Results 

Construction of a minimal model for ERK 
activation by the AT 1A R 

Although some aspects of ERK activation by the AT 1A R are 
understood, the mechanisms controlling the balance and distinct 
kinetic properties of the G protein- and (3-arrestin-dependent 
pathways activating ERK remain largely unknown. To address 
this question, we have developed a deterministic dynamical 
model made of a system of ordinary differential equations 
(ODE). Previous attempts to consider all distinct biochemical 
species and interactions led to combinatorial complexity 
associated with a dramatic increase in the number of unknown 
parameters and, as a consequence, were impractical for 
dynamical simulation (Borisov et al, 2005; Hlavacek et al, 
2006; Birtwistle et al, 2007). Therefore, we chose to construct a 
minimal model encompassing the G protein- and (3-arrestin- 
dependent transduction mechanisms while keeping key mole- 
cules and experimentally measurable variables. 

To build the structure of the model, we established and 
iteratively refined a simplified molecular interaction network 
capturing the available knowledge and experimental data. We 
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generated an initial version of the minimal model amenable to 
numerical simulations and optimization of the unknown 
parameters by training with experimental data. Importantly, 
we found out that with the initial network structure, parameter 
estimation did not lead to simulations satisfyingly fitting 
experimental data. This was not surprising since, as is 
generally the case for biological systems, the knowledge 
available on the signaling network was incomplete and a 
number of competing hypotheses existed. We took advantage 
of this situation by exploring a variety of model structures to 
minimize discrepancies between the simulated model and the 
experimental data. The model structure was reassessed back 
and forth until a satisfactory solution was reached. Impor- 
tantly, the main structural predictions that arose during this 
process were experimentally validated (see below) . The final 
model resulting from this iterative process is a network of 
molecular interactions, which can be formally described either 
using a diagrammatic notation such as Systems Biology 
Graphical Notation (SBGN) (Figure 1; Kitano et al, 2005) or 
equivalently by a set of reaction rules written in the Systems 
Biology Markup Language (SBML) or in BIOCHAM syntax 
(Supplementary Figure SI). 

For simplification purpose, only one cellular compartment 
(volume = 1 pL) was considered. In addition, neither degrada- 
tion nor synthesis of the different molecular species was taken 
into account, and we used mass action laws to model the 
dynamics of all reactions. This led to 11 kinetic reactions, 
depending on 26 kinetic rates (Supplementary Figure S2), and 
7 conservation laws (Supplementary Figure S3) . 

Classical G-protein coupling 

Signal transduction is initiated by hormone (i.e., angiotensin) 
binding to the receptor. At the angiotensin concentration used 
in our experiments (i.e., 100 nM), 99% of the receptor is ligand 
bound. Therefore, for simplification purpose, we aggregated 
the hormone binding process in a single variable called HR. 
The hormone-receptor complex then couples to heterotrimeric 
G protein, ocq subunit is activated (G_a) by exchanging GDP 
with GTP and is subsequently released (Northup et al, 1980; 
Samama et al, 1993; Lefkowitz, 1998). Since this process is 
rapid, transient and dependent upon receptor activation, in the 
model we expressed the whole heterotrimeric G-protein 
activation/deactivation process as G being reversibly activated 
in G_a under the catalytic control of HR. GRK induces 
phosphorylation of the receptor thereby leading to the 
formation of phosphorylated receptor HRP1 (Lefkowitz, 
1998). However, the phosphorylation only partially quenches 
G-protein activation (Benovic et al, 1987). Consequently, 
HRP1 also catalyzes G_a formation in the model. It is also 
well documented that 7TMRs induce G-protein activation 
constitutively (Samama et al, 1993). We took this process 
into account in the model by adding the activation of G 
into G_a in a parallel reaction independent of either HR or 
HRP1 catalysis. 

Desensitization of G-protein activation and recycling of 
the receptor 

Complete desensitization of G-protein activation is achieved 
upon (3-arrestin binding to the GRK-phosphorylated receptor 
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Figure 1 Model for ERK activation by AT 1A R. Red-colored arrow indicates the starting point of the activation process (i.e., ligand binding to the receptor). Green- 
colored species correspond to experimental read-outs. Targeted experimental perturbations are indicated in red. Blue-colored reactions and species correspond to new 
features of the model as compared with classical views in GPCR transduction. Hypothesis made during the modeling process are numbered and circled in blue (1, co- 
existence of two distinct phosphorylated forms of the receptor, HRP1 and HRP2; 2, reversibility of p-arrestin-dependent ERK phosphorylation; 3, p-arrestin-dependent 
ERK phosphorylation undergoes enzymatic amplification; 4, non-phosphorylated ligand-bound receptor still has the ability to recruit p-arrestin 2 and to induce ERK 
phosphorylation; and 5, differential phosphorylation of the ligand-bound receptor by GRK2/3 versus GRK5/6 leading to HRP1 and HRP2). Cell Designer has been used 
to represent the topology of the network in Systems Biology Graphical Notation (SBGN). The following semantic is used: state transition (—□►); proteins ( l ] ): active 
protein (o); simple molecule (O); receptor (^); catalysis (_*_J; association (>-); dissociation HO; phosphorylation (®); Boolean logic gate OR (». Complexes 
are surrounded by a box. 



(Benovic et al 1987; Lefkowitz, 1998). In the model, HRP1 
interacts and forms a complex with either p-arrestin 1 (barrl) 
or p-arrestin 2 (barr2) . Those two complexes are internalized; 
the receptor is then either degraded or dephosphorylated and 
recycled back at the plasma membrane (Lefkowitz and 
Whalen, 2004). As explained above, we did not consider 
synthesis/degradation in the model. The multiple depho- 
sphorylation/recycling steps were shortened in one single 
irreversible reaction departing from each complex and 
providing barrl, barr2 and HR to the system. 

G protein-dependent signaling to ERK 

Once activated, Gaq protein induces a second messenger 
signaling cascade sequentially involving PLC activation, IP3 
and DAG accumulation from PIP2, calcium release and PKC 
activation (Milligan and Kostenis, 2006) . Since we were able to 
measure DAG accumulation and PKC activity, we included 
these variables in the model: all the reactions leading to DAG 
activation/synthesis are modeled as one reaction transforming 
PIP2 into DAG catalyzed by G_a (Figure 1). In turn, DAG 
catalyzes PKC activation in PKC_a. Both DAG accumulation 
and PKC activation have been shown to be deactivated 
through complex mechanisms (Newton, 1995), accordingly 
both reactions were made reversible in the model. Activated 
PKC triggers the activation of many downstream targets 
including Ras and the ERK MAPK signaling module: Raf, 
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MEK, ERK (Wei et al 2003; Luttrell and Gesty-Palmer, 2010). 
ERK phosphorylated through G protein-dependent mechan- 
isms was noted GpERK in the model. Intermediate reaction 
was omitted in the model, and the whole process represented 
as one single reaction: PKC_a catalyzes GpERK formation 
from ERK. This reaction was made reversible to account for 
the activity of dual specificity ERK phosphatases (Dickinson 
and Keyse, 2006). 



p-Arrestin-dependent signaling to ERK 

(3-Arrestin 2 contributes to ERK phosphorylation, indepen- 
dently of G proteins, through the formation of a multiprotein 
complex with the GRK-phosphorylated receptor, Raf and Mek 
(DeFea et al 2000; Luttrell et al 2001; DeWire et al 2007). 
Since this mechanism is known to be differentially affected by 
GRK2/3 or GRK5/6 depletion (Kim et al 2005), we made the 
assumption that the phosphorylated form of the receptor 
engaged into p-arrestin-dependent ERK phosphorylation was 
distinct from the desensitized form. In addition, p-arrestin 2 
depletion has been reported to negatively impact p-arrestin- 
dependent ERK signaling while p-arrestin 1 removal leads to 
increased signaling by this pathway (Ahn et al 2004b). 
Consistently, in the model, a distinct GRK-phosphorylated 
form of the receptor, HRP2, associates with p-arrestin 2 
(barr2) to form the complex HRP2barr2 (Figure 1, hypothesis 
1). HRP2barr2 then promotes the reversible formation of 
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P-arrestin-dependent phosphorylated ERK (bpERK) from the 
cellular pool of ERK (here again Raf and MEK were omitted 
in the model) . GpERK and bpERK were distinguished in the 
model since it is well documented that they exhibit distinct 
kinetics, subcellular localization and downstream targets 
(Ahn et al, 2004a; DeWire et al, 2007; Luttrell and 
Gesty-Palmer, 2010). 

During the parameter optimization process (see below), we 
found that the (3-arrestin-dependent pathway for ERK activa- 
tion needed to be negatively regulated to fit the biological data. 
As shown in Supplementary Figure S4, when rate constants 
kl7 (parameter 39 controlling HRP2 dephosphorylation), k24 
(parameter 46 controlling HRP2barr2 dissociation) or k25 
(parameter 47 controlling bpERK dephosphorylation) are set 
to 0, total pERK formation remains stable after reaching its 
maximum. In the experimental data, however, pERK forma- 
tion clearly decreases between 5 and 90min. As a conse- 
quence, we made the reaction leading to bpERK formation 
reversible (Figure 1, hypothesis 2). In addition, based on the 
actual quantities of receptors (0.08 umol 1 ~ 1 ; Ahn et al, 2004a; 
Kim et al, 2005) of p-arrestin 2 (0.483 umoll" 1 ; Ahn et al, 
2004b) and of phosphorylated ERK (1.265 |imoll" l ; Dupuy 
et al, 2009) measured in HEK293 cells, ERK phosphorylation 
within the p-arrestin scaffold could not work with a 1:1 
stoichiometry. Therefore, we hypothesized the existence of an 
enzymatic amplification within the (3-arrestin scaffold. Using a 
recently validated reverse phase protein array method (Dupuy 
et al, 2009), we were able to experimentally test this 
hypothesis. We measured the molar quantities of phosphory- 
lated MEK and of phosphorylated ERK in the whole cell as well 
as within the (3-arrestin 2 scaffold (i.e., cells treated with a PKC 
inhibitor (Ro-31-8425) to disrupt GpERK) (Supplementary 
Figure S5). In all cases, we found substantially more 
phosphorylated ERK than phosphorylated MEK (5- to 36-fold, 
depending on the conditions). With these data in hand, we 
switched from the equimolar scaffold representation to 
enzymatic catalysis of bpERK formation (Figure 1, hypothesis 
3). Since it is well documented that ERK activated by the (3- 
arrestin-dependent mechanism accumulates within cytosolic 
vesicles, which also contain the receptors and (3-arrestins (Ahn 
et al, 2004a; DeWire et al, 2007), our current data suggest that 
ERK leave the (3-arrestin scaffold once they are phosphorylated 
but remain trapped in the vesicles through a mechanism yet to 
be identified. 



GRK-independent p-arrestin signaling 

Although GRK phosphorylation is well known to significantly 
enhance (3-arrestins' affinity for the receptor (Lefkowitz, 
1998), the main driving force for p-arrestin recruitment and 
activation is the agonist-induced transconformation of the 
receptor itself (Gurevich and Benovic, 1993; Lefkowitz and 
Shenoy, 2005) . Moreover, DeWire et al (2007) have reported 
that AT 1A R A324 with a truncated C-terminal tail is still capable 
of recruiting p-arrestins (albeit weakly) and inducing ERK 
phosphorylation in a p-arrestin-dependent manner. With that 
in mind, we made the assumption that p-arrestin 2 could be 
active even in the absence of GRK phosphorylation. We tested 
this hypothesis by comparing wild-type AT 1A R (WT) with a 
mutant (13 A) in which all the 13 serines and threonines 
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present in the C-tail of the receptor were replaced by alanines 
(Figure 2A) . When expressed in HEK293 cells at similar levels, 
the two receptors behaved differently. In control conditions, 
the 13 A triggered significantly more ERK phosphorylation 
than the WT, which possibly reflects the lack of GRK-mediated 
desensitization (Figure 2B; Supplementary Figure S6). As 
expected, siRNA-mediated p-arrestin 2 knock-down led to a 
dramatic decrease in ERK phosphorylation by the WT, 
especially at late time points (Figure 2C; Supplementary 
Figure S6A and B). Importantly, using the 13 A, p-arrestin 2 
depletion also led to a significant inhibition of the late ERK 
response, strongly supporting our initial assumption that the 
non-phosphorylated receptor could still activate bpERK 
(Figure 2D; Supplementary Figure S6C). As a consequence, 
we added a reversible complexation reaction between HR and 
barr2 to form HRbarr2, which catalyzes bpERK formation 
(Figure 1, hypothesis 4). 



GRK isoforms act differentially at the receptor 

Upon receptor activation, GRK2 and GRK3 are translocated to 
the plasma membrane where they form a complex with the 
free Py subunits of heterotrimeric G protein, making these 
kinases' activation dependent on G protein oc/Py subunit 
dissociation. On the other hand, GRK5 and GRK6 are 
constitutively bound to the plasma membrane and can 
therefore interact with and phosphorylate activated receptors 
independently of heterotrimeric G protein (Lefkowitz, 1998). 
Noteworthy, GRK2/3 and GRK5/6 have been shown to have 
specialized regulatory functions on signaling by various 
7TMRs including the AT 1A R: activation of the p-arrestin 
2-dependent ERK pathway specifically requires GRK5 and 
GRK6 action while second messenger generation is negatively 
regulated by GRK2, but is unaffected by GRK5 or 6. In addition, 
GRK2 and GRK3 negatively impact on ERK phosphorylation 
(Kim etal, 2005; Ren etal, 2005; Gesty-Palmer et al, 2006; Kara 
et al, 2006) . These four GRK isoforms have also been shown to 
participate in receptor phosphorylation, with GRK2 contribut- 
ing the most (Iwata et al, 2005 ; Kim et al, 2005 ; Ren et al, 2005 ; 
Kara et al, 2006). However, it is not known from the literature 
whether GRKs regulate p-arrestin-dependent ERK activation 
by phosphorylating the activated receptor or whether they act 
on an unknown intermediate (Kim et al, 2005; Reiter and 
Lefkowitz, 2006) . This question was central to the construc- 
tion of our model, especially since our experimental data 
showed that p-arrestin signaling also occurs in the absence of 
receptor phosphorylation in the C-tail (Figure 2A-D) . To sort 
out the two possibilities, we compared WT and 13 A AT 1A R 
for their sensitivity to siRNA-mediated GRK depletion. 
Supplementary Figure S7B-D shows representative results 
for siRNA-mediated depletion of the endogenous GRK2, 5 
and 6, respectively. We reasoned that if GRKs acted on an 
intermediate, their depletion would impact ERK activation in 
cells expressing either the 13 A or the WT AT 1A R. On the 
contrary, if they acted by phosphorylating the receptor, then 
their depletion would have no effect on the 13 A. As expected, 
GRK2 and GRK5 differentially regulated ERK response by the 
WT receptor (Figure 2E and F; Supplementary Figure S8A and 
B). In cells expressing the 13 A, however, GRK2 depletion was 
unable to enhance ERK phosphorylation whereas GRK5 
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Figure 2 Role of the AT 1A R C-tail phosphorylation in ERK activation. HEK293 cells were transfected with either WT or a mutant (1 3A) AT 1A R. The indicated siRNAs 
were transfected simultaneously with the receptors. Values in graphs are expressed as percent of control (CTL) stimulated for 5 min with Angll (100 nM). (A) C-terminal 
sequence of the WT (Ser and Thr residues in blue) and 13A (mutated residues in red) AT 1A R. (B) Time course of pERK response by WT versus 13A. (C, D) Effect of 
p-arrestin 2 depletion on pERK response by WT and 13A, respectively. (E, G) Effect of GRK2 depletion on pERK response by WT and 13A respectively. (F, H) Effect of 
GRK5 depletion on pERK response by WT and 13A, respectively. Mean ± s.e.m. from at least four independent experiments with *P<0.05; **P<0.01; ***P< 0.001 
and NS, not significant, as compared with controls. 



knock-down did not modify the ERK response (Figure 2G and 
H; Supplementary Figure S8C and D) . This experiment allowed 
us to better delineate GRKs' action and further supported the 
previously proposed concept of GRK2/3 and GRK5/6 phos- 
phorylating distinct residue combinations on activated 7TMRs 
(Kim et al, 2005; Ren et al, 2005; Gesty-Palmer et al, 2006; Kara 
et al, 2006) (Figure 1, hypothesis 5). Since the available data 
(Kim et al, 2005) did not allow us to discriminate GRK2 from 
GRK3 or GRK5 from GRK6, we aggregated these kinases into 
two variables: GRK23 and GRK56 (Figure 1). In the model, 
HRP1 formation is catalyzed by GRK23 and HRP2 formation 
by GRK56 with HRP1 and HRP2 referring to receptor 
phosphorylated on different serine and/or threonine residue 
combinations. 
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Parameter estimation 

The model comprises a total of 59 parameters whose values 
had to be determined (Table I). Supplementary Figure S9 
shows parameter numbers on the SBGN network representa- 
tion. Among the 59 parameters, 16 (Table I, gray shaded) were 
either taken from the literature (parameters 10, 11, 13, 16, 17, 
18 and 19), set to 0 because known to be strictly agonist- 
induced (parameters 1, 4, 5, 6, 7, 8 and 9) or calibrated using 
DAG and PKC activity FRET assays and preliminary simula- 
tions (parameters 2 and 3). In addition, initial quantities of 
bpERK and GpERK in the four perturbed conditions were 
calculated from the data set (Table I, green-shade, parameters 
49, 50, 52, 53, 55, 56, 58 and 59) (Kim et al, 2005). As the 
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Figure 3 Fitting of the model to the experimental data with the optimized parameter set 4. Experimental data are represented as point ± s.e.m. (in black). All the curves 
were generated through model simulations. Data from (A-D) are from Kim et a/ (2005). Data from (E, F) were generated in living cells using FRET sensors for DAG and 
PKC activity, respectively, after stimulation with 100nM Angll. 



remaining 35 parameters (i.e., kO to k35) were experimentally 
unreachable, we optimized their value by training the model 
with experimental data according to a global optimization 
strategy (Supplementary Figure S10). To limit data hetero- 
geneity, we only considered data collected in a unique cellular 
model (i.e., HEK293 cells). We used kinetic data on three 
variables of the model: (i) previously published angiotensin- 
induced ERK phosphorylation data (Kim et al, 2005), (ii) 
original DAG accumulation and (iii) PKC activity data that we 
measured in real time in AT 1A R expressing HEK293 cells by 
using previously described FRET sensors (Violin et al> 2003). 
DAG and PKC data were available in control conditions (i.e., 
90 min stimulation with 100 nM Angll) whereas ERK data were 
obtained in control as well as in four perturbed conditions (i.e., 
90 min stimulation of cells depleted in (3-arrestin 2, GRK2/3, 
GRK5/6 using siRNA as well as in cells treated with a PKC 
inhibitor. Since they were very close, the ERK data obtained in 
GRK2 and 3 or GRK5 and 6-depleted cells, respectively, were 
averaged. 

The parameter optimization process was carried out in the 
Biochemical Abstract Machine (BIOCHAM) modeling envir- 
onment (Calzone et al, 2006; Fages and Rizk, 2008). This 
software combines Quantifier-Free Linear Logic QFLTL(R) 
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constraints (Rizk et al, 2008) with the Covariance Matrix 
Adaptation Evolution Strategy (CMA-ES), a previously vali- 
dated non-linear optimization technique (Hansen and 
Ostermeier, 2001) to infer parameter values in high dimension. 
Out of 50 optimization runs, we identified 11 sets of parameter 
values that gave model simulations fitting with experimental 
data. Importantly, the initial parameter values used for 
each optimization were independent from the others as they 
were randomly chosen. To further discriminate between the 
parameter sets, we challenged them against published data 
sets that had not been used in the parameter estimation (see 
below). Out of the 11 parameter sets, only the 4 presenting the 
lowest global error (sets 1-4) were able to correctly simulate 
the additional data. Figure 3 shows a fit of the original data set 
with parameter set 4. The simulations were all within the error 
bars of the experimental data except for DAG between 15 and 
40 min for which the simulations were not able to capture the 
observed rebound (Figure 3E). This problem could be due to 
either inappropriate model structure or experimental artifact. 
Further investigations will be necessary to discriminate 
between these two possibilities. Of note, the four sets of 
parameters led to very similar fits to the experimental data set 
(Supplementary Figure Sll). 
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Table I Parameters of the model (set 4) 



Parameter # Description Optimized Value References/assumptions 

Initial quantities 



1 
1 


Activated G protein (G_a) 


u (imoi i 


No constitutive activation at the 
ai i A k ^retrei ana uiauser, zuuyj 


z 


JJ/\Lr 


n HOC! nmnll"^ 

u.uuy |imoi i 


Calibrated according to FRET assays 






and simulations 


a 

D 


/\LLlvc 1 IvLy ^1 J\Ly_d J 


D DD9 nmnl 1 ~~ 1 
U.UUZ |J,lllUl 1 


a 1 1 r\*r^3 tori !irpnrni'nfT tr\ Th'T? Th 1 T tcc^i/c 
LydllUldlcU dLLUlUlllg LU riVLJl doodyb 

and simulations 


4 


HRP1 


0|imoll" 1 


AT 1A R phosphorylation is agonist 
induced (Kim et al, 2005) 


5 


HRPl-p arrestin 1 (HRPl-barrl) 


Oumoll" 1 


P-Arrestin recruitment to AT 1A R is 






0|imoll _1 


agonist induced (Kim et al, 2005) 


6 


HRPl-p arrestin 2 (HRPl-barr2) 


P-Arrestin recruitment to AT 1A R is 
agonist induced (Kim et al, 2005) 


7 


HRP2 


Oiimoll" 1 


ATi A R phosphorylation is agonist 






Oumoll" 1 


induced (Kim et al, 2005) 


8 


HR- p-arrestin 2 (HR-barr2) 


P-Arrestin recruitment to AT 1A R is 






agonist induced (Kim et al, 2005) 


9 


HRP2-p-arrestin 2 (HRP2-barr2) 


0|imoll _1 


P-Arrestin recruitment to AT 1A R is 






0.015 urnoH" 1 


agonist induced (Kim et al, 2005) 


10 


bpERK 


Dupuy et al (2009) 


11 


GpERK 


0.015 limoll" 1 


Dupuy et al (2009) 



Total quantities 








12 


G protein (G + G a) 


k28 


56.99 iimoir 1 




13 


Receptor (HR + HRP1 + HRPl-barrl + 




0.08|imoll" 1 


Ahn et al (2004a,b) 




HRPl-barr2 + HR-barr2 + HRP2 + HRP2-barr2) 








14 


DAG (PIP2 + DAG) 


k29 


1.006 |imoir 1 




15 


PKC (PKC + PKC_a) 


k30 


8.842 iimoir 1 




16 


P-Arrestin 1 (barrl + HRPl-barrl) 




0.858 limoll" 1 


Ahn et al (2004a,b) 


17 


P-Arrestin 2 (barr2 + HRPl-barr2 + 




0.483 (imoll" 1 


Ahn et al (2004a,b) 




HR-barr2 + HRP2-barr2) 








18 


ERK (ERK + bpERK + GpERK) 




4.273 limoll" 1 


Dupuy et al (2009) 


Norms 


19 


ERK norm 




0.013 liamol" 1 


Dupuy et al (2009) 


20 


DAG norm 


k34 


4.12l(imol _1 




21 


PKC_a norm 


k35 


7.211(imol" 1 




Reaction rates 








22 


G auto activation 


kO 


3.11e-4min _1 




23 


Activation of G by HRP1 


kl 


O.OlSliimol^min" 1 




24 


Activation of G by HR 


k2 


7.6l(imol~ 1 min~ 1 




25 


Activation of DAG by G_a 


k3 


4.631|imol~ 1 min~ 1 




26 


Activation of PKC by DAG 


k4 


0.0787 l|imol _1 min " 


i 


27 


Phosphorylation of ERK by PKC 


k5 


2.651|imol~ 1 min~ 1 




28 


Deactivation of G_a 


k6 


5.1 min~ 1 




29 


Deactivation of DAG 


k7 


0.461 min" 1 




30 


Deactivation of PKC 


k8 


1.77 mm" 1 




31 


Dephosphorylation of GpERK 


k9 


3.04 min" 1 




32 


Phosphorylation of HR by GRK23 x GRK23 quantity 


klO 


2.051 min" 1 




33 


Association of HRP1 with p-arrestin 1 


kll 


2.61 1 (imol " 1 min " 1 




34 


Association of HRP1 with p-arrestin 2 


kl2 


2.59 liimol" 1 min" 1 




35 


Dissociation of HRPl-barrl complex 


kl3 


0.0062 min" 1 




36 


Dissociation of HRPl-barr2 complex 


kl4 


0.031 min" 1 




37 


Recycling of HRPl-barrl complex 


kl5 


6.54e-5min" 1 




38 


Recycling of HRPl-barr2 complex 


kl6 


0.0723 mm" 1 




39 


Dephosphorylation of HRP2 


kl7 


0.0665 mm" 1 




40 


Phosphorylation of HR by GRK56 x GRK56 quantity 


kl8 


0.896 mm" 1 




41 


Association of HR with barr2 


kl9 


0.205 liimol" 1 min" 1 




42 


Association of HRP2 with barr2 


k20 


1.041 (imol" 1 min" 1 




43 


Phosphorylation of ERK by HR-barr2 


k21 


4.2e-4 1 (imol " 1 min " 


i 


44 


Phosphorylation of ERK by HRP2-barr2 


k22 


14.44 liimol" 1 min" 1 




45 


Dissociation of HR-barr2 complex 


k23 


1.05 mm" 1 




46 


Dissociation of HRP2-barr2 complex 


k24 


0.347 mm" 1 




47 


Dephosphorylation of bpERK 


k25 


0.762 min" 1 




Inhibited PRC condition 








48 


PKC total quantity (PKC + PKC_a) 


k31 


0.0018 limoll" 1 




49 


Initial bpERK quantity 




0.0037 ^imoll" 1 


Kim et al (2005) 


50 


Initial GpERK quantity 




0.0037 ^imoll" 1 


Kim et al (2005) 
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Table I (Continued) 



Parameter # Description Optimized Value 


References/ assumptions 


Inhibited p-arrestin 2 condition 

51 p-arrestin 2 total quantity k32 


1.12e-4 umoll 1 




52 Initial bpERK quantity 

53 Initial GpERK quantity 


8.2e-4 umoll" 1 
8.2e-4 umoll" 1 


Kim et al (2005) 

JYLlll ei LLL ^ZUUj J 


Inhibited GRK23 condition 

54 Phosphorylation rate of HR by GRK23 x GRK23 quantity k26 


1.087 mm" 1 




55 Initial bpERK quantity 

56 Initial GpERK quantity 


0.0182 umoll" 1 
0.0182 umoll" 1 


Kim et al (2005) 
Kim et al (2005) 


Inhibited GRK56 condition 

57 Phosphorylation rate of HR by GRK56 x GRK56 quantity k27 


6.13e-4 min" 1 




58 Initial bpERK quantity 

59 Initial GpERK quantity 


0.0055 umoll" 1 
0.0055 umoll" 1 


Kim et al (2005) 
Kim et al (2005) 





Gray-shaded parameters were either taken from the literature (parameters 10, 11, 13, 16, 17, 18 and 19), set to 0 as they are strictly agonist-induced (parameters 1, 4, 5, 
6, 7, 8 and 9) or calibrated using FRET assays and simulations (parameters 2 and 3). Green-shaded parameters were calculated from our experimental measurements. 
All the remaining parameters (k0-k35) were computationally optimized. All the concentrations were calculated with the assumption that cellular volume was 1 pL. 
Initial bound receptor quantity (HR) was calculated from AT 1A R expression levels of 200-300 fmol per mg of protein measured in HEK293 cells in our transfection 
conditions (Ahn et al, 2004a,b). The fraction of occupied receptor was calculated for different hormone concentrations using mass action law with a Kd of 1.6 nM. 



Parameters identifiability and sensitivity 

We next compared the optimized values of each parameter 
across the four best parameter sets and determined their 
respective sensitivities. To address this question, the value of 
each individual parameter was scanned across 14 logs (from 
10 ~ 9 to 10 5 ) and the impact on the global error (a function of 
the differences between simulated and observed values) was 
determined. When plotted against each other, these values 
provided a clear picture of each parameter's identifiablity (as 
estimated by the convergence of the four independent 
optimization runs), and sensitivity/robustness. Interestingly, 
the parameter values showed a clear convergence (i.e., 
maximum 1 log difference across the four parameter sets) for 
20 out of 35 of the optimized parameters (k2, k4, k5, k6, k8, k9, 
klO, kl2, kl7, kl8, kl9, k20, k22, k24, k25, k26, k28, k29, k30 
and k34) when the four sets were compared (Figure 4; 
Supplementary Figure S12; Supplementary Table SI). In 
addition, two parameters (k3 and kl6) were convergent for 
three parameter sets, the fourth one being tolerant (i.e., 
parameter values that do not significantly impact on the global 
behavior of the model) . Five parameters (k7, kll , kl 3 , kl4 and 
k35) were convergent for three parameter sets. Finally, eight 
parameters (kO, kl, kl5, k21, k23, k27, k31 and k32), even 
though not converging to a similar value, actually presented 
large and overlapping ranges of tolerance (Figure 4) . Together, 
these data strongly suggest that our global optimization 
strategy led to the efficient identification of most parameters 
within the network. Parameter set 4 was used for figures but 
the four sets give very similar behavior in all envisaged 
situations. 



Validation against independent experimental 
data sets 

Using the same experimental design, Kim et al (2005) reported 
that GRK2 or 3 overexpression led to decreased ERK 
phosphorylation upon angiotensin stimulation. Remarkably, 

8 Molecular Systems Biology 2012 



the observed pattern was reminiscent of the one obtained in 
GRK5- or 6-depleted cells. Oppositely, GRK5 or 6 overexpres- 
sion caused an increase in phosphorylated ERK at all time 
points, similarly to what was observed in GRK2- or 3 -depleted 
cells. When increases in GRK23 (Figure 5A) or GRK56 
(Figure 5B) quantities were simulated, the model predicted 
phosphorylated ERK patterns which matched those reported in 
Figure 5 of Kim et al (2005). 

Next, we took advantage of another interesting set of data 
generated in the same experimental conditions as those we 
have used (Ahn et al, 2004a) . Using either a mutant version of 
AT 1A R (DRY-AAY) or a biased ligand (SII), these authors 
showed that, in the absence of G-protein coupling, ERK was 
still triggered and displayed kinetics virtually identical to the 
ones obtained in cells treated with a PKC inhibitor. We 
therefore inhibited G proteins in the model (Figure 5C) and 
that resulted in simulated levels of phosphorylated ERK similar 
to PKC inhibition shown in Figure 3B and in perfect agreement 
with the data reported by Ahn et al (2004a) in their Figure 2. 

Finally, another paper reported that, in the case of AT 1A R, 
P-arrestin 1 acted as a dominant negative of p-arrestin 2 for 
ERK phosphorylation upon angiotensin stimulation (Ahn et al, 
2004b). By simulating the effects of either p-arrestin 1 or 
P-arrestin 2 depletion on phosphorylated ERK response as a 
function of angiotensin concentration (Figure 5D), we also 
observed that a potentiating effect was triggered by p-arrestin 1 
depletion albeit the amplitude of the effect was markedly 
weaker than the one observed experimentally in Figure 1C of Ahn 
et al (2004b) . The opposite effect of p-arrestin 1 and 2 was more 
pronounced in simulations performed at 5 min with parameter 
set 3 (Supplementary Figure S13A) as well as at 10 and 30 min 
with parameter set 4 (Supplementary Figure S13C and D). It 
should be noted however that, in their experiments, Ahn et al 
analyzed ERK phosphorylation from cytosol-enriched extracts, 
predominantly containing p-arrestin-regulated phosphorylated 
ERK, rather than from whole cell ly sates. Since the model 
simulates ERK phosphorylated in the whole cell, this might 
account for the observed difference in the kinetics. 
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Figure 4 Comparison of four optimized parameter sets. The optimized value of each parameter is figured by a dark color line and compared side-by-side for the four 
parameter sets. The tolerance ranges (i.e., all the values of one parameter that induce less than a three-fold increase in the global error of the model) are figured by color- 
shaded bars. Set 1 is blue, set 2 is yellow, set 3 is red and set 4 is green. Parameters are numbered as presented in Table I. 



Interestingly, our simulations provide quantitative estima- 
tions of the degrees of the different inhibition, depletion or 
overexpression tested (Table I; Supplementary Table SI). In 
particular, the values found for the perturbed parameters were 
dramatically diminished ( > 99 % inhibition) in all cases except 
GRK23 for which the reduction was only about 50% of the 
control value (Table I, compare parameters 32 and 54). This 
result suggests that, unlike GRK5 and GRK6, GRK2 and GRK3 
exert additive non-redundant effects on (3-arrestin-induced 
ERK in vivo. We further assessed the model validity by doing a 
series of simulations using selected perturbations (e.g., 
Supplementary Figures S14, SI 5, SI 6 and SI 7). We found that 
the model generally behaved as expected, predicting activa- 
tion/deactivation half-lives consistent with independent 
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experimental measurements reported in the literature 
(Supplementary Table S2; Ahn et al, 2004a; Rajagopal et al, 
2006; Violin etal, 2006; Lohse et al, 2008; Vilardaga, 2010; Poll 
et al, 2011). This further validates the model and the global 
optimization strategy we have used for its parameterization. 



Pivotal role of GRK23 in the control of ERK 
activation by p-arrestin 2 

Interestingly, we also made the unexpected finding that, in the 
absence of GRK23, the amount of HRP2 and HRP2barr2 nearly 
doubled, suggesting that bpERK might be substantially 
increased upon GRK23 depletion (Supplementary Figure 
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Time (min) Log (Angll) 



Figure 5 Model validation against published independent data sets. Control conditions are presented in blue; gradual increases or decreases (as indicated) of the 
targeted parameter are represented using different gradations of purple. (A) Simulation of GRK23 overexpression (GRK23 quantity multiplied by 2, 5 or 10). 
(B) Simulation of GRK56 overexpression (GRK56 quantity multiplied by 2, 5 and 1 0). (C) Inhibition of G-protein activation (total quantity of G protein divided by 2, 5 or set 
to 0). (D) Simulation of either (3-arrestin 1 (red curve) or p-arrestin 2 (green curve) depletion on pERK response at 5 min as a function of angiotensin concentration. 



S17A versus D). We then determined that when GRK23 was 
depleted, the model predicted an increase in both second 
messenger (DAG) (Figure 6A; Supplementary Figure SI 5) and 
total pERK (Figure 6C) compared with control conditions. The 
classical paradigm would have explained the increased ERK 
response as resulting from a lack of G-protein desensitization 
in the absence of GRK23. However, total pERK was minimally 
affected by PRC blockade when combined with GRR23 
depletion. Simulations clearly predicted that the contribution 
of GpERR is decreased instead and that bpERR activation is 
massively amplified (Supplementary Figure S18). Importantly, 
these predictions were corroborated by experimental data. To 
reflect second messenger response, inositol uptake was 
measured in the presence of increasing concentrations of 
angiotensin. As in the simulation, second messenger accumu- 
lation was higher in GRR2-depleted cells than in control 
(Figure 6B) . Moreover, we observed that ERR phosphorylation 
was strongly increased and that the PRC inhibitor had 
very limited effect in the absence of GRR2 (Figure 6D; 
Supplementary Figure SI 9). Although we were not able to 
achieve double GRR2 and GRR3 siRNA-mediated depletion 
experimentally, the experimental observations matched the 
model simulations with good accuracy. 

To explore the possibility that GRR23 exerts a strong 
negative effect on the p-arrestin-dependent pathway through 
its competition with GRR56 for receptor phosphorylation, we 
then simulated the combined effect of GRR56 depletion and 
PRC inhibition on ERR phosphorylation. At short stimulation 
time, GRR56 depletion had very little effect on the amount of 
ERR phosphorylation as compared with control conditions. 
However, PRC inhibition combined with GRR56 depletion led 
to complete abolition of ERR response whereas significant, 
albeit reduced, ERR response was predicted by the model in 
PRC inhibited control siRNA-transfected cells (Figure 6E). 
Again, these model simulations were closely corroborated by 
the experiments. We observed that ERR phosphorylation upon 
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exposure to increasing doses of angiotensin was similar in 
control and GRR5- or GRR6-depleted cells. Moreover, the PRC 
inhibitor completely blocked ERR phosphorylation in GRR5- 
or GRR6-depleted cells while leading to only partial inhibition 
in control cells (Figure 6F and G; Supplementary Figure S20). 
These data clearly indicate that GRR5 or 6 depletion directs 
ERR signaling exclusively toward the G protein-dependent 
transduction mechanism. Interestingly, siRNA-mediated 
depletion of either GRR5 or 6 is sufficient to completely 
redirect the ERR signal. This suggests that the two kinases may 
have to phosphorylate the receptor on distinct residues (or 
combination of residues) to impart (3-arrestin-mediated 
signaling competence to this receptor. 

These data put forward the intriguing hypothesis that GRR2 
might exert a strong negative effect on p-arrestin-dependent 
ERR activation. This suggests that this kinase, rather than 
being restricted to the desensitization of G-protein signaling as 
classically thought, might exert an equally important dampen- 
ing action on p-arrestin signaling through its competition with 
GRR5 and 6 for receptor phosphorylation. To determine 
whether or not this mechanism exists in an authentic cell 
context, we used rat primary vascular smooth muscle cells 
(VSMCs) as a model. It had been previously shown that, in rat 
VSMC, ERR is phosphorylated upon angiotensin exposure 
through a bimodal mechanism involving PRC and p-arrestins 
similar to the one modeled in HER293 cells (Ahn et dl, 2009; 
Rim et al, 2009). Our data clearly demonstrate that, as 
observed in HER293 cells, activation of endogenously 
expressed AT 1A R in GRR2-depleted VSMC leads to a sub- 
stantial increase in PRC-independent pERR (Figure 7A-C). 

In addition, we provide evidence that a similar mechanism 
also operates with other 7TMRs. Indeed, the endogenously 
expressed P2AR, as well as overexpressed vasopressin V2 
(V2R), FSHR and neurokinin 1 (NR1R) receptors all displayed, 
albeit to different extent, increased PRA or PRC-independent 
pERR in response to GRR2 depletion in HER293 cells 

© 2012 EMBO and Macmillan Publishers Limited 



Competing GRKs balance 7TMR signaling 
D Heitzler et al 



Simulated 



A 

A 






200 i 


1% 




'E 

13 


150 - 






£5 






100 - 


!5 




c5 




O 


50 - 


< 




Q 






0 - 




-1 


c 






200 -I 


<z> 




'E 




13 


150 - 


>, 








100 - 








50 - 






LU 






0 - 




-1 


E 






100 i 






'E 

13 


80- 








60- 








40- 






cr 


20- 


LU 




Q. 






0- 




-1 



- CTL siRNA 




Log (Angll) 

— CTL siRNA 

— GRK2/3 siRNA 

— GRK2/3 siRNA + PKC inhibitor 




Log (Angll) 

- CTL siRNA 

- GRK5/6 siRNA 

- PKC inhibitor 

GRK5/6 siRNA + PKC inhibitor 




Log (Angll) 



B 



Experimental 



CO x 

o £ 



oc E 

LU 3 



CO 



LU 



* I 

cc e 

LU x 

1 2 



CC 
LU 



250 -i 


■ CTL siRNA 








a GRK2 siRNA 






200 - 






150 - 








100 - 




■ 


50 - 








0 < 





-11 -10 -9 -8 -7 
Log (Angll) 



■ CTL siRNA 
a GRK2 siRNA 

• GRK2 siRNA + PKC inhibitor 



-6 




125 
100 



-10 -9 -8 
Log (Angll) 

CTL siRNA 
GRK5 siRNA 

CTL siRNA + PKC inhibitor 
GRK5 siRNA + PKC inhibitor 




Log (Angll) 



125 
100 
75 
50 
25 
0 



■ CTL siRNA 
^ GRK6 siRNA 

GRK6 siRNA + PKC inhibitor 
' CTL siRNA + PKC inhibitor 




-10 -9 -8 
Log (Angll) 



Figure 6 Role of GRK isoforms in regulating second messenger generation and subsequent ERK signaling by AT 1A R. Simulations used the optimized parameter set 
with 0.1 Limoll" 1 as the initial GRK23 quantity in GRK23-depleted condition. (A) DAG response at 10 min as a function of angiotensin in GRK2/3-depleted (purple) 
compared with control (blue). (B) Cells transfected with AT 1A R and either control (blue) or GRK2 (purple) siRNAs. After 10 min of exposure to increasing angiotensin 
doses, inositol phosphate levels were measured. (C) Simulated pERK response at 5 min as a function of angiotensin concentration in control (blue), GRK2/3-depleted 
(purple), or GRK2/3-depleted and PKC-inhibited by Ro-31-8425 (black) conditions. (D) Cells were transfected with AT 1A R and either control (blue) or GRK2 (purple and 
black) siRNAs. Serum-starved cells were pre-incubated with DMSO (blue and purple) or with Ro-31-8425 (Ro, black) and stimulated for 2 min with increasing doses of 
angiotensin. (E) Simulated pERK response at 2 min as a function of angiotensin concentration in either control (blue), GRK5/6-depleted (purple), PKC-inhibited (red) or 
GRK5/6-depleted and PKC-inhibited (green) conditions. (F, G) Cells were transfected with AT 1A R and control (blue and red), GRK5 (F) or GRK6 (G) (purple and green) 
siRNAs. Serum-starved cells were pre-incubated with DMSO (blue and purple) or Ro-31-8425 (black and green) and stimulated for 2 min with increasing doses of 
angiotensin. Signals were normalized according to the amount of total ERK. Mean ± s.e.m. from at least five independent experiments. 



(Figure 7D-G; Supplementary Figure S21) . Together, these data 
strongly support the notion that the GRK-regulated balance 
between G protein- and p-arrestin-dependent signaling that is 
highlighted here, is not restricted to AT 1A R overexpressed in 
HEK293 cells. On the contrary, the uncovered mechanism 
likely applies to other 7TMRs and to physiologically relevant 
situations. 



Discussion 

Seven transmembrane receptor research is challenged with the 
question of how G protein- and p-arrestin-dependent transduc- 
tion mechanisms coordinate their actions to activate common 
downstream signaling intermediates with distinct temporal 
and spatial patterns. In this study, we combined computational 
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Figure 7 Evidence for GRK2-mediated regulation of the balance between G protein and p-arrestin through endogenously expressed AT 1A R in primary VSMCs and 
through four other GPCRs in HEK293 cells. (A) Depletion of endogenous GRK2 in control versus GRK2 siRNA-transfected VSMC. (B) Representative western blots of 
primary rat VSMC transiently transfected with either control (CTL) or GRK2 siRNAs. Serum-starved VSMC was pre-incubated with DMSO or with a PKC inhibitor 
(G66983 10liM) and stimulated for the indicated time with angiotensin (1 ul/l). (C) Quantification of ERK phosphorylation kinetics in control versus GRK2-depleted 
VSMC pretreated with either DMSO or a PKC inhibitor. (D) Time course of ERK phosphorylation upon stimulation of endogenously expressed P2AR with isoproterenol 
(10u.M) in HEK293 cells transfected with CTL or GRK2 siRNA and pre-treated with DMSO or with a PKA inhibitor (H-89, 10 jaM). (E) Time course of ERK 
phosphorylation upon stimulation of transiently transfected V2R with dvd-AVP (0.1 |iM) in HEK293 cells co-transfected with CTL or GRK2 siRNA and pre-treated with 
DMSO or with a PKA inhibitor (H-89, 1 0 uM). (F) Time course of ERK phosphorylation upon stimulation of transiently transfected FSHR with FSH (3 nM) in HEK293 cells 
co-transfected with CTL or GRK2 siRNA and pre-treated with DMSO or with a PKA inhibitor (H-89, 10 liM). (G) Time course of ERK phosphorylation upon stimulation of 
transiently transfected NK1 R with Substance P (0.1 liM) in cells co-transfected with CTL or GRK2 siRNA and pre-treated with DMSO or with a PKC inhibitor (Ro-31- 
8525, 1 u.M). In all experiments, signals were normalized according to the amount of total ERK. Results are expressed as mean ± s.e.m. from at least three independent 
experiments. 



modeling and experiments to address this question using the 
well-documented dual ERK phosphorylation by AT 1A R as a 
paradigm. We brought together previously published experi- 
mental data with our own new data and we assembled an 
original model inference and parameter optimization strategy 
to develop a dynamical model. We chose to build an abstracted 
model rather than a detailed mechanistic model encompassing 
all the molecular details because it had already been 
established that such an approach leads to combinatorial 
increase of complexity which cannot currently be computa- 
tionally handled (Borisov et al, 2005; Hlavacek et al, 2006; 
Birtwistle et al, 2007). This model construction and analysis 
led to new insights into the regulation of ERK by the AT 1A R. 

A common limitation encountered when modeling intracel- 
lular signaling networks is that the number of measurable 
read-outs is very limited with respect to the number of 
molecular intermediates actually involved in the pathways to 
be modeled (Swameye et al, 2003 ; Raue et al, 2009) . This raises 
problems of identifiability since many of the parameters have 
to be determined through an optimization process dependent 



on the available experimental data (Raue et al, 2009). 
Therefore, a limited number of read-outs implies a relative 
lack of constraint on the parameter optimization process, 
which in turn precludes the convergence of parameter sets 
separately optimized. In the present study, since we were 
limited to three read-outs (i.e., pERK, DAG and PKC), we had to 
approach the model identifiability problem in an innovative 
manner. We achieved a global fitting strategy that encom- 
passed both control conditions and a number of targeted 
inhibitions within the network. We found out that, using this 
approach, the constraints on the parameter search were 
significantly higher. In addition to curve fitting, we took 
advantage of the capabilities of BIOCHAM to express in Linear 
Time Logic QFLTL(R) to combine a number of qualitative 
constraints further decreasing the degrees of freedom. Using 
this approach, starting from 50 independent sets of randomly 
chosen parameter values, we obtained four independent 
parameter sets that fitted the experimental data with very 
low error. Interestingly, the four parameter sets presented a 
high degree of convergence with most of the 35 optimized 
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parameters being identifiable. The global optimization strat- 
egy validated in the present study allowed us to parameterize a 
model comprising 18 molecular species by actually measuring 
only three of them over time. In its principle, this global 
approach is versatile and could be adapted to many of the 
dynamical modeling needs that are encountered in systems 
biology. 

As anticipated, having the ability to simulate a dynamical 
model proved to be very efficient to decipher our signaling 
system and its key regulatory features. However, we also 
realized that the model-building process itself was very 
insightful and led to important findings on the network 
structure. More specifically, we found out that the goodness of 
fit resulting from the parameter optimization process was a 
very efficient way to rank different network structures and to 
ultimately infer a model from prior knowledge and data. This 
is important since model inference remains a major challenge 
in systems biology (Nelander et al, 2008; Xu et al, 2010). 

Besides AT 1A R, a growing number of 7TMRs have been 
reported to activate ERK phosphorylation through dual G 
protein/ p-arrestin-mediated transduction mechanisms, which 
share many attributes with the AT 1A R- dependent signaling 
network modeled in this study (Kim et al, 2005; Gesty-Palmer 
et al, 2006; Kara et al, 2006; Shenoy et al, 2006; Luttrell and 
Gesty-Palmer, 2010). One aspect that seems to vary from one 
receptor to another is the role of (3-arrestin 1. Some receptors 
require both (3-arrestin isoforms to activate ERK while for 
others, including AT 1A R, (3-arrestin 1 antagonizes the action of 
(3-arrestin 2 (Ahn et al, 2004a) . Despite of this, we observed 
that one novel mechanism predicted by the model (i.e., GRK2- 
mediated inhibition of (3-arrestin-dependent signaling) seems 
to apply to multiple 7TMRs (i.e., AT 1A R, (32AR, V2R, FSHR and 
NK1R) and different cellular context (i.e., HEK293 versus 
primary VSMC) . This supports the notion that our core model 
has a generic value, as it could predict a signaling mechanism 
used by other 7TMRs and/or in other cellular contexts. Of 
course, to achieve optimal results, the parameters would 
have to be optimized, on a case by case basis, using specific 
data sets. 

The main feature revealed by our modeling approach is that 
ERK activation by AT 1A R is tightly controlled by the 
antagonistic actions of the two GRK subfamilies, namely 
GRK2 and 3 on the one hand and GRK5 and 6 on the other. The 
model made the assumption that distinct phosphorylated 
forms of the receptor (i.e., HRP1 and HRP2) were formed as a 
consequence of GRK2/3- or GRK5/6-mediated actions. This 
assumption was supported by the indirect experimental 
evidence presented in Figure 2 as well as by reports from the 
literature. For instance, this idea is consistent with previous 
results showing that the presence or absence of serine and 
threonine clusters in the receptor C terminus regulates the 
affinity of (3-arrestin recruitment and the pattern of intracel- 
lular trafficking for a wide number of 7TMRs (Oakley et al, 
2000, 2001). It is also in agreement with the fact that different 
GRK subtypes differently regulate G protein- and (3-arrestin- 
dependent signaling by various 7TMRs (Iwata et al, 2005; Kim 
et al, 2005; Ren et al, 2005; Kara et al, 2006; Shenoy et al, 
2006). As the present study was carried out, three papers were 
published that strongly support the prediction of our model. A 
study revealed that the dynamically regulated site-specific 
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phosphorylation of the CXCR4 by multiple kinases leads to 
both positive and negative regulation of CXCR4 signaling 
(Busillo et al, 2010). A separate study demonstrated that the 
M3 -muscarinic receptor C terminus is preferentially phos- 
phorylated on specific sites when stimulated with full versus 
partial agonists (Butcher et al, 2011) . Finally, the detailed GRK- 
dependent phosphorylation 'bar code' occurring at the (32AR 
has been deciphered in another independent study (Nobles 
et al, 2011). Interestingly, six GRK2- and two GRK6-dependent 
non-overlapping sites were identified. Both GRK2- and GRK6- 
dependent sites contributed to desensitization and internaliza- 
tion. However, only GRK6-dependent sites were required for 
(3-arrestin-mediated ERK activation. Remarkably, carvedilol, 
a previously characterized weak (3-arrestin-biased ligand 
(Wisler et al, 2007), only induced phosphorylation on the 
GRK6 sites, whereas isoproterenol, a full agonist at the (32AR, 
triggered phosphorylation on both GRK2 and GRK6 sites. 
Importantly, phosphorylation on the GRK6 sites was increased 
in GRK2-depleted cells, providing direct experimental evi- 
dence for the existence of a competition between the two GRKs 
for receptor phosphorylation. 

In addition, our simulations and experimental data suggest 
that these GRKs operate as the major control gate regulating 
ERK activation in time and space. In the absence of GRK5 or 6, 
there is virtually no (3-arrestin-dependent ERK activation 
whereas, more surprisingly, in the absence of GRK2, 
(3-arrestin-dependent ERK activation clearly overwhelms G 
protein-mediated ERK. In the absence of GRK action, both 
pathways co-exist but the potential for regulation is much 
reduced. This reveals a new and unexpected function of GRK2: 
beyond its well-documented role in G-protein desensitization, 
this kinase strongly dampens (3-arrestin signaling. Eventually, 
it is the interplay between GRK2 and GRK5/6 that defines the 
duration and intensity of each transduction mechanism. We 
propose that GRKs endow ERK activation (and potentially 
other signaling pathways) by 7TMRs with adaptability. The 
balance in the action of both GRK subfamilies will vary with 
the receptor (i.e., different phosphorylation patterns present in 
the receptors' C terminus), the cell type and the cellular 
context (i.e., differences in stoichiometries and subcellular 
localizations) . 

The fact that GRK2 on the one hand and GRK5 and 6 on the 
other are both required for the regulatory process to operate 
emphasizes the importance of tightly controlling the balance 
between the G protein- and the (3-arrestin-dependent transduc- 
tion mechanisms. Importantly, this is consistent with the 
mounting evidence that they both lead to distinct biological 
outcomes (DeWire et al, 2007; Luttrell and Gesty-Palmer, 
2010). 

Importantly, (3-arrestin-dependent signaling is not limited to 
ERK activation. Instead, (3-arrestins are now considered to be G 
protein-independent signal transducers acting as multifunc- 
tional scaffolds that interact with many protein partners 
(Xiao et al, 2007) and protein kinases, thereby leading to 
the phosphorylation of numerous intracellular targets 
(Christensen et al, 2010; Xiao et al, 2010) . The model developed 
in the present study provides a robust core encompassing the 
proximal transduction events involving G proteins and 
(3-arrestins as well as their regulatory mechanisms. We believe 
that our model could provide a solid core module that could 
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accommodate the incorporation of additional G protein- and/ 
or p-arrestin-dependent signaling cascades. 

The concept of 'biased agonism' that has recently emerged 
in pharmacology defines another important application field 
for our model. In particular, one class of 'biased ligands' that 
has the ability to alter the balance between G protein- 
dependent and p-arrestin-dependent signal transduction 
engenders a lot of interest both in academia and in 
pharmaceutical industry (Reiter et al, 2012). In that context, 
a real added value of our model is that it provides a clear and 
rational mechanistic basis to this novel class of compounds. In 
that, the model proposed here represents a first step toward 
rational discovery of G protein and p-arrestin-biased drugs 
active at different 7TMRs. In particular, the model will help 
identify better suited methods allowing the screening and/or 
characterization of new biased compounds. 

Materials and methods 

Materials 

Radiolabeled [ 125 I] Tyr 4 -angiotensin II (Angll) was obtained from 
Perkin-Elmer Life Sciences (Boston, MA). PRC (Ro-31-8425 and 
G66983) and PRA (H89) inhibitors were purchased from Calbiochem 
(Darmstadt, Germany) . GeneSilencer reagent was from Gene Therapy 
Systems (San Diego, CA). All other reagents were purchased from 
Sigma (St. Louis, MO). Expression plasmids encoding the rat 
hemagglutinin epitope-tagged (HA-) AT 1A R (Wei et al, 2003), human 
FLAG-V2R (Ren et al, 2005), rat FLAG-FSH-R (Reiter et al, 2001) and 
human FLAG-NR1R have been used. All the 13 Ser and Thr residues 
present in the C-tail of the rat AT 1A R have been mutated to alanine 
(13 A) using the QuickChange site-directed mutagenesis kit (Strata- 
gene, La Jolla, CA). 



Synthesis of siRNAs 

Chemically synthesized double-stranded siRNAs with 19 nt duplex 
RNA and 2 nt 3' dTdT overhangs were purchased from Dharmacon 
(Lafayette, CO) or Xeragon (Germantown, MD) in deprotected and 
desalted form. The siRNA sequences used in the present study to target 
the different human p-arrestins and GRKs have been previously 
validated (Ahn et al, 2003; Kim et al, 2005; Ren et al, 2005; Kara et al, 
2006): (3-arrestin 1, 5'-AGCCUUCUGCGCGGAGAAU-3' (positions 
441-459 relative to the start codon); (3-arrestin 2, 5'-GGACCGCAAAGU 
GUUUGUG-3' (positions 150-168); GRK2, 5' -GAAGUACGAGAAGCUG 
GAG-3' (positions 270-288); GRK 5, 5'-GCCGUGCAAAGAACUCUUU-3' 
(positions 408-426); GRK6, S'-CAGUAGGUUUGUAGUGAGC-S' (posi- 
tions 726-744). The siRNA sequence used to target rat GRK2 in VSMC 
cells was S'-AAGAAAUAUGAGAAGCUGGAG-S' (positions 270-288). A 
non-silencing RNA duplex (S'-UUCUCCGAACGUGUCACGU-S') was 
used as a control. 



Cell culture and RNA transfection 

HEK293 cells were cultured in minimum Eagle's medium supplemen- 
ted with 10% fetal bovine serum and 1% penicillin/streptomycin. 
HEK293 cells were transiently transfected as described (Kim et al, 
2005; Kara et al, 2006). Briefly, 30-40% confluent, slow growing early 
passage ( < 10) cells in 100 mm dishes were transfected simultaneously 
with 20 |ig of siRNA and plasmid encoding the AT 1A R (WT: 0.5 (ig; 13A: 
2.0 ng), V2R (1.0 ng), FSHR (1.0 ng) or NK1R (1.0 ug) using the 
GeneSilencer Transfection reagent (Gene Therapy System, San Diego, 
CA) . Forty-eight hours after transfection, cells were divided into poly- 
D-lysine-coated 12- well plates (Becton Dickinson Labware, Bedford, 
MA) for receptor binding, or 12- well plates to prepare cell extracts. 
AT 1A R expression was determined by radioligand binding assays, as 
described previously (Laporte et al, 1996), and was 300-600 fmol per 



mg of protein in all experiments. Expression levels of V2R, FSHR and 
NK1R in CTL versus GRK2 siRNA- transfected cells were monitored by 
ELISA using anti-FLAG antibody. In the GRK siRNA experiments, the 
amount of plasmid transfected was adjusted to obtain equivalent 
receptor levels at the plasma membrane as in the controls. 

Rat VSMCs were prepared from aorta of male Sprague-Dawley rats 
by enzymatic digestion and maintained in DMEM supplemented with 
1 % glutamine, 10 % fetal bovine serum and 1 % penicillin-streptomy- 
cin as previously reported (Kim et al, 2009). In all, 80-90% confluent, 
slow growing early passage (<6) cells in 100 mm dishes were 
transfected with 20|ig of either control or rat GRK2 siRNA using 
60 ul of Lipofectamine 2000 transfection reagent (Invitrogen, Carlsbad, 
CA) . Forty-eight hours after transfection, cells were starved with media 
with 0.1% bovine serum albumin for least 24 h before stimulation 
(Ahn et al, 2009; Kim et al, 2009). 



Preparation of cell extracts and immunoblotting 

Three days after transfection, cells were starved for at least 6h in 
serum-free medium before stimulation. After stimulation, cells were 
solubilized in 2 x SDS-sample buffer (pH 6.8), followed by sonication 
or boiling. Equivalent amounts of proteins were separated by SDS/ 
PAGE on Tris • glycine polyacrylamide gels (Invitrogen) , transferred 
onto nitrocellulose, and immunoblotted with rabbit polyclonal 
antibodies. The phospho (1:2000 dilution) and total (1:6000 dilution) 
ERK antibodies were from Cell Signaling Technology (Beverly, MA) 
and Upstate (Charlottesville, VA), respectively. Endogenous p-arrest- 
ins and GRKs were detected as described (Ahn et al, 2003; Kim et al, 
2005). Chemiluminescent detection was performed using the Super- 
Signal West Pico or West Femto reagent (Pierce, Rockford, IL) and was 
quantified by densitometry with a Fluor-S Multilmager (Bio-Rad, 
Hercules, CA) or with ImageMaster ID Elite version 4 software 
(Amersham Biosciences, Pittsburgh, PA) . 



FRET 

HEK293 cells expressing the ATi A R (600 fmol per mg of protein) were 
transiently transfected with either DAG (DAGR) or PKC (CKAR) 
plasmid-encoded FRET sensors (Violin et al, 2003). Twenty-four hours 
after transfection, cells were split and seeded on imaging dishes (BD 
Biocoat, BD Biosciences, San Jose, CA, USA) . Twenty-four hours later, 
cells were serum starved for 90 min and imaged in the dark on a 37°C 
temperature-controlled stage. Cells were imaged on a Leica DM IRBE 
(Leica Microsystems, Wetzlar, Germany) microscope with a CoolSnap 
fx cooled charge-coupled device camera (Roper Scientific, Ottobrunn, 
Germany) controlled by METAFLUOR 7.5 (Universal Imaging Corpora- 
tion, Downingtown, PA). Dual emission ratio imaging used a 436DF10 
excitation filter, a 436-510 DBDR dichroic mirror, and 480AF30 and 
550AF30 emission filters for CFP and YFP, respectively. The filters were 
alternated by a Lambda 10-2 filter changer (Sutter Instruments, 
Novato, CA). Exposure time was 100-500 ms, and images were taken 
every 15 s. Typically, 15 sensor-positive cells and 5 non-specific areas 
were chosen in the field of the microscope. Baseline signals at the 
wavelength corresponding to YFP and CFP were recorded for 5 min 
before the hormone was added to the dish. Signals were recorded for 
the next 90 min. Fluorescent image background corrections were 
carried out by subtracting the intensity of non-specific areas. The FRET 
ratios were calculated for each individual cell. The mean values ± " 
s.e.m. corresponding to individual cells from three different plates 
were plotted. 



RPPA 

Serial dilutions of phosphorylated human recombinant ERK2 (Pro- 
teinkinase, Biaffin GmbH & Co KG, Germany) and of phosphorylated 
human recombinant MEK1 (Millipore, Billerica, MA) were prepared in 
RPPA lysis buffer. We previously determined that >99% of the 
purified pERK2 was indeed dually phosphorylated (Dupuy et al, 2009) . 
In the case of pMEKl, based on the provider certificate of analysis, we 
made the assumption that 86.5% of the preparation was dually 
phosphorylated MEK1. For cell lysates (1.7 mgrnl - 1 concentration in 
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average) , three two-fold serial dilutions were done in RPPA lysis buffer. 
Samples were printed on nitrocellulose-coated slides (Fast Slides, 
Whatman, Maidstone, UK) using a 32-pin manual arrayer (Glass Slide 
Micro arr ay er, VP478, V&P Scientific, San Diego, CA), immunodetected 
and revealed in near infrared according to previously reported 
procedures (Dupuy etal, 2009). Anti-pERK (Cell Signaling Technology 
Inc., Beverly, MA; 1/400), anti-pMEK (Cell Signaling Technology Inc.; 
1/300) and anti-ERK2 (Santa Cruz Biotechnology Inc., Santa Cruz, CA; 
1/1000) antibodies were used. Finally, slides were air dried by 
centrifugation at 2500r.p.m. for 25min and scanned at 700 nm with 
an Odyssey IR imaging system (LI-COR Biosciences, Lincoln, NE) at a 
42-um resolution. Scanned images of arrays were analyzed with 
GenePix Pro 6.0 software (Molecular Devices, Sunnyvale, CA) as 
described (Dupuy et al, 2009). 



Inositol phosphate determination 

HEK293 cells expressing the AT 1A R were transiently transfected in 
100 mm dishes with siRNA targeted against the different GRK 
subtypes. Two days post transfection, cells were plated onto poly-D- 
lysine-coated 12-well plates (BD Biosciences Labware, San Jose, CA). 
To assay for inositol phosphate production, cells were incubated 
overnight at 37°C in labeling medium (2.5 uCi of myo-[ 3 H] inositol in 
0.5 ml of minimum Eagle's medium with 10% fetal bovine serum/ 
well). Cells were washed with Hank's balanced salt solution contain- 
ing 20 mM LiCl for 15min at 37°C and then treated with agonist in 
Hank's balanced salt solution for lOmin. The reactions were 
terminated by the addition of perchloric acid, and total inositol 
phosphates were isolated by anion-exchange chromatography on 
Dowex AG1-X8 columns as described previously (Cotecchia et al, 
1992). 



Statistical analysis 

Statistical analysis of the data was performed with two-way analysis 
of variance to compare entire curves (GraphPad Software Inc., 
San Diego, CA). 



Parameter optimization 

For parameter estimation, we used a newly established parameter 
optimization approach combining quantifier-free linear time logic 
QFLTL(R) (Fages and Rizk, 2008; Rizk et al, 2009) for expressing 
qualitative constraints, with the covariance matrix adaptation evolu- 
tion strategy CMA-ES (Hansen and Ostermeier, 2001) for non-linear 
optimization. The whole parameter optimization process was carried 
out in the BIOCHAM modeling environment (Calzone et al, 2006; 
Fages and Rizk, 2008). Fifty global optimizations were made each 
starting from a random set of parameters and therefore independent 
from each other. Eleven optimizations reached global errors <0.1. 
When challenged against the additional published data (Figure 5), 
only the four parameter sets presenting the lowest global errors led to 
simulations that matched with the experimental data. These four 
parameter sets were selected for further analysis (Supplementary Table 
SI). The ODE model and best parameter set are available from the 
Biomodels database under the ID: MODEL1012080000 (http://www. 
ebi.ac.uk/biomodels/) . 



Additional constraints to the objective function 

To infer the unknown (kinetic and quantitative) parameter values from 
the experimental data obtained under various conditions, we 
expressed the objective function to be minimized in temporal logic 
properties with quantifier free LTL(R) formulae (QFLTL(R); Rizk et al, 
2008) . This formalism allowed us to define, in addition to curve fitting, 
a number of constraints that were integrated to the objective function: 
(i) each parameter value had to be within 0 and 200; (ii) HRP2 
recruited at least five times more barr2 than HR (parameter 42 > 5 x 
parameter 41); (iii) HRP1 recruited more barrl than barr2 (parameter 
33 > parameter 34); (iv) inhibited parameter < control parameter 



(siRNA depletions of GRK23, GRK56 and barr2, PRC inhibition); 
(v) constant boundaries (as defined by standard deviations) imposed 
for DAG (between 5 and 35min) and PRC (between 20 and 80min) 
values instead of all the discrete time points; and (vi) basal levels 
of phosphorylated ERK were low (if HR = 0, then l%<GpERK + 
bpERK<6%). 



Sensitivity analysis 

We carried out a sensitivity analysis to test the robustness and the 
identifiability of each of the computationally optimized parameter 
across the four best parameter sets. The value of each individual 
parameter was scanned across 14 logs (from 10 ~ 9 to 10 s ) and the 
impact on the global error of the fitting to the learning data set was 
determined. Variations of the global errors were then plotted as a 
function of each parameter's value. To assess the convergence of the 
four independent optimization runs, the errors were normalized and 
compared on the same plot for each parameter. In addition, the relative 
'tolerance range' was calculated for each parameter as being the range 
of a parameter value that increases the global error less than three-fold. 
These data were then compared for all the optimized parameters 
across the four parameter sets using box plots. Box plots present for 
each parameter the optimized value within its 'tolerance range'. 



Supplementary information 

Supplementary information is available at the Molecular Systems 
Biology website (www.nature.com/msb). 
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